#### Convergence checks (z, t, checks take some time)
library(coda)

# do Geweke and Heidel-Welch diags
load("post.uds.RData")

# total number of params in model
nparams <- ncol(post.uds$z[[1]]) + ncol(post.uds$gamma[[1]]) +
ncol(post.uds$sigmasq[[1]]) + ncol(post.uds$t[["arat"]][[1]]) +
ncol(post.uds$t[["blm"]][[1]]) + ncol(post.uds$t[["bollen"]][[1]]) +
ncol(post.uds$t[["freedomhouse"]][[1]]) +
ncol(post.uds$t[["hadenius"]][[1]]) +
ncol(post.uds$t[["pacl"]][[1]]) + ncol(post.uds$t[["polity"]][[1]]) +
ncol(post.uds$t[["polyarchy"]][[1]]) + ncol(post.uds$t[["prc"]][[1]]) +
ncol(post.uds$t[["vanhanen"]][[1]])

# heidelberger and welch
hs <- heidel.diag(post.uds$sigmasq)
hg <- heidel.diag(post.uds$gamma)
hz <- heidel.diag(post.uds$z)
harat <- heidel.diag(post.uds$t[["arat"]])
hblm <- heidel.diag(post.uds$t[["blm"]])
hbollen <- heidel.diag(post.uds$t[["bollen"]])
hfhouse <- heidel.diag(post.uds$t[["freedomhouse"]])
hhadenius <- heidel.diag(post.uds$t[["hadenius"]])
hpacl <- heidel.diag(post.uds$t[["pacl"]])
hpolity <- heidel.diag(post.uds$t[["polity"]])
hpolyarchy <- heidel.diag(post.uds$t[["polyarchy"]])
hprc <- heidel.diag(post.uds$t[["prc"]])
hvanhanen <- heidel.diag(post.uds$t[["vanhanen"]])
# 12-rater model
#hmunck <- heidel.diag(post.uds$t[["munck"]])
#hmainwaring <- heidel.diag(post.uds$t[["mainwaring"]])

h <- c(hs[,1], hg[,1], hz[,1], harat[,1], hblm[,1], hbollen[,1],
hfhouse[,1], hhadenius[,1], hpacl[,1], hpolity[,1], hpolyarchy[,1],
hprc[,1], hvanhanen[,1])
# 12-rater model
#h <- c(h, hmunck[,1], hmainwaring[,1])

print(sum(h==0))

# geweke
gs <- geweke.diag(post.uds$sigmasq)
gg <- geweke.diag(post.uds$gamma)
gz <- geweke.diag(post.uds$z)
garat <- geweke.diag(post.uds$t[["arat"]])
gblm <- geweke.diag(post.uds$t[["blm"]])
gbollen <- geweke.diag(post.uds$t[["bollen"]])
gfhouse <- geweke.diag(post.uds$t[["freedomhouse"]])
ghadenius <- geweke.diag(post.uds$t[["hadenius"]])
gpacl <- geweke.diag(post.uds$t[["pacl"]])
gpolity <- geweke.diag(post.uds$t[["polity"]])
gpolyarchy <- geweke.diag(post.uds$t[["polyarchy"]])
gprc <- geweke.diag(post.uds$t[["prc"]])
gvanhanen <- geweke.diag(post.uds$t[["vanhanen"]])
# 12-rater model
#gmunck <- geweke.diag(post.uds$t[["munck"]])
#gmainwaring <- geweke.diag(post.uds$t[["mainwaring"]])

g <- c(gs$z, gg$z, gz$z, garat$z, gblm$z, gbollen$z, gfhouse$z,
ghadenius$z, gpacl$z,
  gpolity$z, gpolyarchy$z, gprc$z, gvanhanen$z)
# 12-rater model
#g <- c(g, gmunck$z, gmainwaring$z)

print(sum(abs(g)>2))

save(g, h, gs, gg, gz, garat, gblm, gbollen, gfhouse, ghadenius,
gpacl, gpolity, gpolyarchy, gprc, gvanhanen, hs, hg, hz, harat, hblm,
hbollen, hfhouse, ghadenius, hpacl, hpolity, hpolyarchy, hprc,
hvanhanen, file="geweke-heidel.RData")

# clean up
rm(list=ls())
